TD1 - Statistique Spatiale

library(tidyverse)
library(ggforce)
library(gstat)
library(sf)


dta_full_sf <- st_read(dsn = "swiss_full/swiss_rain_full.shp") 
dta_sf <- st_read(dsn = "swiss/swiss_rain.shp")
bords_sf <- st_read(dsn = "swiss/bords.shp") 
elev_sf <- st_read(dsn = "swiss/elev.shp") 

dta_sf |> st_coordinates() |> bind_cols(dta_sf) |> ggplot() + 
    geom_sf(data = bords_sf, linewidth = 0.1) +
    geom_circle(aes(x0=X, y0=Y, r=pluies/100, fill = pluies)) +
    scale_fill_viridis_c(name = "Pluies", option = "mako", direction = -1) 


dta_full_sf |> st_coordinates() |> bind_cols(dta_full_sf) |> 
    ggplot() + 
    geom_sf(data = bords_sf, linewidth = 0.1) +
    geom_circle(aes(x0=X, y0=Y, r=pluies/100, fill = pluies)) +
    scale_fill_viridis_c(name = "Pluies", option = "mako", direction = -1) 

elev_sf |> ggplot() + 
    geom_sf(aes(col = elevation)) + 
    geom_sf(data = bords_sf, fill = NA, linewidth = 1.5) + 
    theme_minimal() +
    theme(legend.position = "right") +
    scale_color_viridis_c(name = "Z", option = "mako", transform= "log10")

nuees_dta <- do.call(rbind, lapply(1:nrow(dta_sf), function(i) {
    xy <- st_coordinates(dta_sf)
    x0 <- xy[i, 1]
    y0 <- xy[i, 2]
    pluie0 <- dta_sf$pluies[i]
    dta_sf |> rownames_to_column('Id')  |> 
        mutate(delta_x = xy[,1] - x0, delta_y = xy[,2] - y0, ecart = pluies - pluie0) |> 
        filter(Id > i) |> 
        mutate(h = sqrt(delta_x^2 + delta_y^2), ecart2 = ecart^2 / 2) |> 
        select(Id, h, ecart2)
}))

nuees_dta |> ggplot() + geom_point(aes(x=h, y=ecart2)) + xlim(c(0, 120000))

vario.cloud <- variogram(pluies~1, data = dta_sf, cloud = TRUE)
vario.cloud |> ggplot(aes(x=dist, y=gamma)) + geom_point() + ggtitle("Nuée variographique")

vario.iso <- variogram(pluies~1, data = dta_sf)
v.fit <- fit.variogram(vario.iso, vgm(model = "Exp", 15000, 60))
vario_fit_exp <- variogramLine(v.fit, maxdist = 120)

ggplot(data = vario.iso, aes(x=dist, y=gamma)) + 
    geom_point() + geom_line(data= vario_fit_exp , aes(x=dist, y=gamma))

Kfull <- krige(formula = pluies ~ 1, locations = dta_sf, newdata = dta_full_sf, model = v.fit)
Kfull |> ggplot() + geom_sf(aes(fill = var1.pred)) +
    geom_sf(data = dta_sf, aes(col = pluies)) +
    scale_fill_viridis_c(option = "mako", transform = "log10")